subroutine define_parameter
use var
use type
implicit none
real(dp) :: dr, NP_volume
integer :: i,j,k,label

  tnodes = dime(1)*dime(2)*dime(3)

  df = dble(box)/dble(dime-1) ! distance between two nodes
  dr = df(1) ! Correct if NP density equal in all direction

  write(*,"(A,3f7.3)") " Distance between two nodes = " ,df
  write(222,"(A,3f7.3)") " Distance between two nodes = " ,df
  
  df = 1/df ! Invers of df
  
  df2 = df*0.5 ! Invers of 2*df
  
  ddf = df*df

  nodes = 0
  bulk = .false.
  do k=1,dime(3)
     do j=1,dime(2)
        do i=1,dime(1)

           label = 1 + (i-1) + (j-1)*dime(1) + (k-1)*dime(1)*dime(2)
           if (number_NP>0) then
           if (channel) then
              if (boundary(i,j,k)==.false..and.boundary_NP(i,j,k)==.false..and.NP(i,j,k)==.false.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           elseif (droplet.or.cylinder) then
              if (boundary(i,j,k)==.false..and.boundary_NP(i,j,k)==.false..and.NP(i,j,k)==.false..and.drop(i,j,k)==.true.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           endif
           else
           if (channel) then
              if (boundary(i,j,k)==.false.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           elseif (droplet.or.cylinder) then
              if (boundary(i,j,k)==.false..and.drop(i,j,k)==.true.) then
                 bulk(i,j,k) = .true.
                 nodes = nodes + 1
              endif
           endif
           endif
        enddo
     enddo
  enddo
  if (droplet) deallocate(drop)
  if (number_NP>0) deallocate(NP)

  ! differential 
 ! nodes = tnodes - bnodes - sum(NPnodes) - sum(NPbnodes)
  NP_volume = number_NP*((4.0/3.0)*Pi*((radius)*dr)**3)
  if (channel) then
     Volume = (box(1)*box(2)*box(3)) - NP_volume ! For cube
     ad0 = 2*(box(1)*box(2)) / bnodes
  elseif (droplet) then
     Volume = (4.d0/3.d0)*Pi*((center(1)-2)*dr)**3 - NP_volume / 2.0d0
     ad0 = (4.d0*Pi*((center(1)-2)*dr)**2 -  number_NP*Pi*((radius)*dr)**2) / bnodes
  elseif (cylinder) then
     Volume = (Pi*((center(1)-2)*dr)**2)*box(3)
     ad0 = (2.d0*Pi*((center(1)-2)*dr))*box(3) / bnodes
  endif
  vd = Volume / nodes
  Print *, vd, ad0
! Tylor   
!  vd = box(1)*box(2)*box(3)/((dime(1)-2)*(dime(2)-2)*(dime(3)-2))
  if (number_NP>0) then
     if (channel)   ad1 = number_NP*(Pi*((radius)*dr)**2) / sum(NPnodes)
     if (droplet)   ad1 = number_NP*(Pi*((radius)*dr)**2) / sum(NPnodes) / 2.0d0
     if (cylinder)  ad1 = number_NP*(2*Pi*((radius)*dr)*box(3)) / sum(NPnodes)
  else
  ad1 = 0
  endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Eneargy Scale

Escale = Elastic_L1 * nematic_l * 1.0d-9 ! (J/m)*(nm)*(10^-9)


 print *, "first Elastic constat (J/m)", Elastic_L1
 print *, "Nematic coherance length (nm)", nematic_l
 print *, "A                       (J/m3)",  Elastic_L1/(nematic_l * 1.0d-9)**2
 Print *, "Surface Energy          (J/m2)",gama(1)*Elastic_L1/(nematic_l * 1.0d-9),gama(2)*Elastic_L1/(nematic_l * 1.0d-9)
 Print *, "Escale                        ", Escale
 print *, "  "
 print *, "Total nodes         = ", tnodes 
 print *, 'Slab boundary nodes = ', bnodes
 print *, 'NP boundary nodes   = ', sum(NPbnodes)
 print *, 'Nodes inside NP     = ', sum(NPnodes)
 print *, 'MC move nodes       = ', nodes
 write(*,'(40A)') ' #################################################################################################'
 write(*,'(40A)') '    CPU   STEP         GL           ELASTIC         SURFACE0          SURFACE1          TOTAL'
 write(*,'(40A)') ' -------------------------------------------------------------------------------------------------'
 
 write(222,*) "first Elastic constat (J/m)", Elastic_L1
 write(222,*) "Nematic coherance length (nm)", nematic_l
 write(222,*) "A                       (J/m3)",  Elastic_L1/(nematic_l * 1.0d-9)**2
 write(222,*) "Surface Energy          (J/m2)",gama(1)*Elastic_L1/(nematic_l * 1.0d-9),gama(2)*Elastic_L1/(nematic_l * 1.0d-9)


 write(222,*) "Total nodes         = ", tnodes 
 write(222,*) 'Slab boundary nodes = ', bnodes
 write(222,*) 'NP boundary nodes   = ', sum(NPbnodes)
 write(222,*) 'Nodes inside NP     = ', sum(NPnodes)
 write(222,*) 'MC move nodes       = ', nodes
 write(222,'(40A)') ' #################################################################################################'
 write(222,'(40A)') '    CPU   STEP         GL           ELASTIC         SURFACE0          SURFACE1          TOTAL'
 write(222,'(40A)') ' -------------------------------------------------------------------------------------------------'
end subroutine  define_parameter


